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NOMENCLATURE 


ADC 

N 

2N + 1 

f 

f 


f 

s 



s 


G(f) 

Af 

h(n), h , H 
n n 

h L (n) 

h A (n) 

V n > 

h R ( n ) 

D 


Analog -to-digital converter 
Integer associated with number of weights 
Number of weights used to filter 
Frequency variable, usually hertz (Hz) 

Cutoff frequency 
Termination frequency 

Sampling frequency, number of samples per second 
Time between samples 

Gain variable associated with frequencies in describing a 
filter 

Difference between termination frequency and cutoff 
frequency 

Filter weights 

Subscript refers to weight type; the example is a low- 
pass filter 

All-pass filter weights 
Band -pass filter weights 
High-pass filter weights 
Band-rejection filter weights 
Digital data 
Derivative weights 
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HR 
PR 
R, R 1 
F, F 1 


NOMENCLATURE (Concluded) 

Real component of response 
Imaginary component of response 
Normalized frequency variable; that is, f/f 

s 

Frequency variable for plotting in cycles per second 
(cps) 
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TECHNICAL MEMORANDUM X-64697 

DDP-116 GENERAL DIGITAL FILTERING 
INTRODUCTION 

The purpose of this program is to offer a general purpose digital 
filtering capability on the Statistical Analyzer (STAN) system (Appendix A) . 
The filtering methods used in this program are those methods developed by 
Martin-Graham, as discussed in Reference 1. The following filtering capa- 
bilities are available. 


1. All Pass 

3 

2. Low Pass 

3. High Pass 

4. BandPass 

5. Band Rejection 

6. Derivative 


INPUT 


The program possesses capabilities to filter either raw, unconverted 
digital data or engineering unit data. Input tapes can assume any of the 
following forms: 

1. Astrodata ana log -to -digital converter (ADC) 

2. Scientific Data Systems 

3. MSFC format 

4. S&E-COMP-RRF binary tape (same as item 3 except header 
record is missing) 


5. Cycle counter 



The Astrodata ADC, the cycle counter, and the Scientific Data Systems 
produce raw data that can be converted to engineering units using either sine 
or step calibrations. Descriptions of the above tape formats may be obtained 
from the authors. 

CALIBRATION 


The maximum (C 2 ) and minimum (C t ) calibration count levels are 
computed for the sine calibration by 



1.414 



( 1 ) 



and 



N 

£ x i 


Ft 1 


(3) 


where 


X. ■=' the data count in a calibration file 

i 

i = 1, 2, . . . , N , 


and 


N = the number of data points in a calibration file 


The slope is computed by 


S = 



(4a) 


where E 2 and Ej are maximum and minimum engineering units* When a 
sine wave is used for calibration, 



(4b) 


and the offset is computed by 


0 = E 2 - S(C 2 + C) 


(5) 


The raw data points (X^) are converted to engineering units (d.) by 


d. = S • X. + 0 

l i 


( 6 ) 


The maximum (C 2 ) and minimum (Cj) calibration count levels are computed 
for the step calibration. First the maximum count level (C 2 ) is considered. 
The mean (c) and standard deviation (C ) count level for the maximum 
calibration level are computed by 


i N 

c " ■ h 2 c i 

i=l 


i = 1, 2, ... , N, N « 200 


(7) 


•tJFi 


i = 1, 2, ... , N, N « 200 


( 8 ) 
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c = c + c , 

max s 


( 9 ) 


and 


C . = C - C 

mm s 


( 10 ) 


All C. are interrograted for the value (C 2: C. & C . ), and any C. 
i & ' max i min' J i 

outside these limits is rejected. All data within these limits are accepted 

and equations (?) through (10) are repeated until C si. When this condi- 

S’ 

tion is met, the minimum calibration level Cj is computed identically to C 2 . 

After C 2 and Cj are computed, equations (4) through (6) are used to 

compute d.. 

i 

FILTERING 
Filter Weights 

Generally the number of filter weights (2N + l) is chosen such that a 
given amplitude error is maintained over the intervals 0 s f < f and 

f -s f < f /2, where f is the frequency variable, f is the cutoff frequency, 
t s c 

f is the termination frequency, and f is the sampling frequency (Fig. l) . 
l s 



Figure 1. Filter definition. 
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Error Analysis . The error is defined as the difference between the 
ideal filter curve and the actual response curve in the intervals 0^ f^ f 

c 

and f ^ f^. The errors discussed refer to the maximum errors and occur at 

no more than two points. The average error and the expected error are much 
less than the maximum errors discussed. The percent of maximum error is 
defined in terms of the following parameters: 


N — where 2N + 1 is the number of filter weights. 

At — interval between consecutive data points (independent variable) . 

Af — filter termination frequency, f , minus filter cutoff frequency, 

f . 1 

c 


Figure 2 illustrates these results. The table inset in Figure 2 is used to 
compute N, the number of weights being 2N + 1, as follows: 

1. The percent of maximum amplitude error as a function of NAtAf 
is shown in Figure 2. This curve will be used by the analyst to determine the 
optimum selection of N and Af. 

2. Normally At is a fixed parameter; i. e. , the sampling rate should 
be at least four times the bandwidth of the unfiltered function. Using 40 000 
samples per second as the maximum sampling rate, the following values will 
be realized: 


2.5 x 10 5 < At< 1 
f = 10 000 Hz , 


and 

f = 20 000 Hz 

v 
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PERCENT MAXIMUM ERROR 



NAtAf 


Figure 2. Error curve. ' ( ■- 

3. Af is chosen such that all unwanted data are rejected. 

4. N is chosen such that an acceptable error is maintained. The 
operator must judiciously select N to minimize computing time and still 
maintain an acceptable error. 
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5. The optimum error value'for most applications will be 1.5 per- 
cent for which 


0.95 

N = (round N up to the nearest integer) , (ll) 

where 

At = — - — , 

sampling rate 

A f =f t .£ c , 
f = termination frequency , 

and 

f c = cutoff frequency 

All-Pass Filter Weights 

•V 

All -pass filter weights h(n) yield 


h(n) 


0; -N ^ n :S -1 
1 ; n = 0 
0; 1 ^ n < N 


( 12 ) 


Low-Pass Filter Weights 

A test for a singular point occurrence is performed: 


__ 

= TAtRAf) * 


(13) 
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If n is an integer, then 


(14) 


(15) 


There are 2N + 1 weights (h n> h , , h^, ... , h^ ). Equation 

(15) is evaluated for all values of n= 1, 2, 3, ... , N, noting that h(n) = 
h(-n). This yields 2N weights and h(o) = At (f + f) , thus producing the 
central value of the 2N + 1 values. C 

The set of weights h (n) constitutes the weights for a low -pass filter 

J-i 

with the following characteristics: 


h(nAt) = At 


Af 7 7T \ 


For all nonintegers, 


h ( n) = - — 
L v ' 2irn 


sin (27rf nAt) + sin ( 27rf nAt) 
X c 

1-4 (Af) 2 (nAt) 2 


G(f) 


1; 0 < f < f 

c 


J 

'(f -y' 

. 4 

2 J COS 

n 

ti— i 

<1 

J 



0i£ t s fs ( f s - V 

High -Pass Filter Weights 


The following equations illustrate high-pass filter weights: 

h^( 0) = 1 - h^O) for n = 0 
h H (n) = -h L (n) for n ± 0 



(16) 


(17) 
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When combining two sets of weights by equation (17), an equal number of 
weights must be considered by both filters. 

Band-Pass Filter Weights 

The band-pass filter weights h^(n) are defined as the corresponding 
differences between two low-pass weight sets; that is, 

h B (n) = h L (n) - h^n) . (18) 

This has the same effect as geometrically subtracting the filter gain (Fig. 3). 

Band -Rejection Filter Weights 

The band-rejection filter weights ^ BR ( n ) are computed by negating 

all the band-pass weights except the central weight, which is subtracted from 
unity: 



Figure 3. Band -pass filter. 
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y»)-.h A (n) - bgW , 
or 

V*> ‘ - h a (n! • 

except 


V (0) = 1 - h B (0) • 

APPLICATION OF WEIGHTS 

After the filter type is selected and the weights are computed, 
normalized for the low -pass filter only: 

h L (n) 

V = -N ‘ 

X h L (n) 

df -N 

Next, the weights H (n) are applied to the linearized data, with K 

K 

L, BR designating the filter type, 

N 

%+l = H K ( n ) d [(N+n)+l] ’ - 

ITF; ~N 

where 

d (N+n+l) = dl » d2 ’ •" ’ d (n+l) ’ * ’ d (2N+l) * 
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(19a) 


(19b) 


they are 


( 20 ) 


A, B, 


( 21 ) 



The smoothing process is initiated by assuming data values to the left equal 

to those on the right. For example, if raw data starts at d , then d = 

d T _ , where K= 1, , N. s Sf x 

s-K • 

A similar procedure is used for the end points d ; d __ = d 
. „ „ e e-K e+K 


DERIVATIVE WEIGHTS 


If derivative weights are desired the following procedure applies. 
For the central weight. 


D(0) = 0 


( 22 ) 


and, with the exception of a possible singular weight, all other weights are 


D = h 
n n 


[1 - 12(AfAtn) 2 ] - [f At cos(2iTf Atn) + f At cos(2Trf Atn) ] 

t X c c 

nAtll - 4( AfAtn) 2 ] 


(23) 


where 

n = 1, 2, . . . , N 


The negative weights become 


D = -D (24) 

-n n 

/ ' 


The singular derivative weight is computed when 


n 


0.5 

(At) (Af) 


(25) 
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omputed, 


the filter configuration is 


( 27 ) 


( 28 ) 


( 29 ) 


( 30 ) 



TRANSFER FUNCTIONS 


Options are available to plot the selected filter transfer functions or 
response curves. The phase and gain are computed and plotted as a function 
of frequency: 


HR = 


2N+1 / s 

Yj h. cos\[i - (N + l) ] 2 ttR 1 > 

i=l 1 ' ’ 


(31) 


and 


PR = 


2N+1 / 

/. h. sin<[i - 
i=l 1 * 


(N + 1)] 27TR 1 


) ■ 


(32) 


The gain G(f) is computed; 


G(f) = (1 + R) 2 + (PR) 2 


(33) 


The phase response <p is computed; 


<P = tan 1 (PR/HR) 


(34) 


The gain and phase are computed at the discrete points; 


R = 0, 0.01, 0.02, ... , 0.51 


(35) 


and are plotted at the discrete frequency points 


F = Of , O.Olf , 0.02 f , ... , 0. 51 f 
s s s s 


(36) 
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The abscissa is scaled by computing R 1 and F 1 ; 


R 1 = RK + AtC 


( 37 ) 


and 

F 1 = FK + C , F = n(O.Ol) and n = 0, 1, ... , 51 , (38) 


where K is a multiplication factor and C is a constant to shift the response 
curve from zero frequency. 

For example, assume that f = 100, then the response curve will be 

s 

plotted at the discrete points 0, 1, 2, ... , 51; however, if better resolution 
is desired, equations (37) and (38) apply. Assuming the area of interest is 
over the interval 


5S f< 10 


let 


C = 5 


and 

K = 1 . 

Then, 

F l = (0 + 5) , (0.1 + 5) , (0.2 + 5) (5.1 + 5) 


The ordinate for the phase is plotted over the interval -180° ^ ^ 180° . 
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DATA RESULTS 


The plots on the following pages illustrate the various filter configura- 
tions as applied to a square wave that is a composite of odd harmonics extending 
to the upper frequencies limited by the Nyquist frequency. 

CONCLUSION 


The DDP-116 General Digital Filtering program is a general purpose 
filter program designed to run under the STAN operating system. The program 
is designed to interface with all elements of the Data Reduction Branch to pro- 
vide general purpose filtering services. 
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ID SQUARE ME AS. ID. El 1-1 CAL. RANGE —100.000 TO 100.000 

NO. WEIGHTS -80 SING. POINT - 40 CUTOFF FREQUENCY - 100.000 TERMINAL FREQUENCY - 200.000 
ALL-PASS FILTER BANDWIDTH - 200.000 SAMPLE RATE - 3000/mc 



< 

0 

D 

01 

Z 



TIME (nc) 



TIME (nc) 


All -Pass Filter Output and Input 
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ID SQUARE MEA8.ID. E11-1 CAL. RANGE *>-100.000 TO 10&000 

NO. WEIGHTS" 30 SING. POINT - 40 CUTOFF FREQUENCY «* 100.000 TERMINAL FREQUENCY - 209.000 
BAND-PASS FILTER BANDWIDTH - 280.080 SAMPLE RATE - 8000/we 










INPUT DATA FILTERED DATA 









High-Pass Filter Input and Output 








3SVHd 4 NIVO 


NO. WEIGHTS - 80 SING. POINT - 40 CUTOFF FREQUENCY - 100.000 TERMINAL FREQUENCY • 200.000 
LOW-PASS FILTER DERIVATIVE BANDWIDTH - 200.000 SAMPLE RATE - 8000/sac 

RESPONSE CURVES FACTOR - 8000.000 CONSTANT - 0.000 



0 1000 2000 3000 4000 

FREQUENCY (Hi) ' 


Response Curve for a Low-Pass Filter Derivative 
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APPENDIX ; 

STATISTICAL ANALYZER SYSTEM 

The Statistical Analyzer system can be used as a general purpose 
computer system and can also perform the following special operations: 

1. Auto-Correlation 

2. Cross -Correlation 

3. Power Spectra (Harmed) 

4. Cospectra and Quadspectra (Harmed) 

5. RMS Amplitude versus Frequency 

6. Transfer Function Gain and Phase 

7. Coherence Function 

. 8. Single and Joint Probability Distributions 

i 

9. Mean, Variance, Skewness, and Kurtosis 
10. Statistical Confidence Tests 

Functional Description (Fig. A-l) 

The Statistical Analyzer system is comprised of two computational 
devices, a DDP-116 Honeywell computer and an Astrodata convolver; four 
input/output and memory devices; three Control Data 607 tape drives; a 
Datadisc disc memory; a Control Data 405 card reader; a Video Systems 
display unit; and the logic required to interface these devices. 

The computer transmits instructions and data to the devices via the 
interface logic. When a device is ready to participate in a data transfer 
with the computer, the interface logic transmits a data request to the computer. 
In response to this request, one or more data characters are transferred 
between the computer and the requesting device via the interface logic. 
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Figure A-l. Statistical Analyzer system, overall block diagram. 




For each data transfer, an address in the computer memory must be 
specified. Blocks of addresses are specified in terms of starting and ending 
addresses to be associated with transfers between the computer and each of 
the other devices. Having received the proper addresses, the interface logic 
stores and processes these address data in order to specify addresses for 
individual transfers in consecutive order until the block of assigned addresses 
has been used. 

The interface provides various control signals to the devices as 
required to implement instructions received from the computer, synchronizes 
data transfers, formats data as necessary, and summarizes status signals 
received from the devices to simplify monitoring by the computer. 

The special operations performed by the system are accomplished 
through use of the Astrodata convolver. The convolver is an extremely high 
speed arithmetic device and special purpose computer with its own core 
memory. It can respond to the following instructions: 

1. Load 

2. Unload 

3. Auto -Correlate 

4. Cross -Correlate 

5. Compute Power Spectral Density 

6. Multiply and Sum 
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